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Computation of the temperature dependence of the heat capacity of complex molecular systems 
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We propose a new method for computing the temperature dependence of the heat capacity in complex molec- 
ular systems. The proposed scheme is based on the use of the Langevin equation with low frequency color 
noise. We obtain the temperature dependence of the correlation time of random noises, which enables to model 
the partial thermalization of high-frequency vibrations, which is a pure quantum effect. By applying the method 
to carbon nanotubes, we show that the consideration of the color noise in the Langevin equation allows to 
reproduce the temperature evolution of the specific heat with good accuracy. 
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I. INTRODUCTION 



The pronounced temperature dependence of the specific 
heat c{T) in molecular systems is a pure quantum effect. It 
is well known that in the absence of any critical point, the in- 
crease of the temperature is accompanied by a smooth rise of 
the specific heat c{T), while by decreasing the temperature, 
c(r) tends to zero. The first explanation for this phenomenon 
was given by Einstein a century ago HJ ■ This quantum effect 
results from the fact that at low temperatures, high frequency 
vibrations are partly frozen while low frequency vibrations 
are fully excited. It is of course impossible to explain the ef- 
fect of partial thermalization within the framework of classical 
physics. In fact, the use of the usual Langevin equation with 
white noise leads to a uniform thermalization of all modes and 
the specific heat is practically insensitive to the temperature 
(in classical physics, any temperature dependence of the ther- 



mal capacity results from non-linearity effects). On the other 
hand, it is possible to mimic the partial thermalization effect 
in the framework of the Langevin description if one consid- 
ers, instead of a white noise, a low frequency color noise with 
a temperature dependent frequency spectrum. To this end, it 
is enough to introduce a random noise with a finite correlation 
time ic > during which the noise "remembers" its last real- 
ization. In this study, we obtain the temperature dependence 
of this correlation time, which allows a coiTect modeling of 
the partial thermalization of vibrations in the system. 

The paper is organized as follows. We consider in Sec. I a 
harmonic oscillator and investigate its dynamics described by 
the Langevin equation with color noise. We obtain the tem- 
perature dependence of the correlation time of random noises, 
which enables us to efficiently model the partial thermaliza- 
tion of high-frequency vibrations. We next evaluate in Sec. 
II the effect of non-linearities on the accuracy of our com- 



putational method and propose a general scheme to compute 
the heat capacity for many-body systems using a color noise. 
The final section is devoted to the application of the proposed 
scheme to a Hamiltonian carbon nanotube model. It is shown 
in this section that the temperature evolution of the specific 
heat computed with the Langevin equation with color noise 
agrees well with the result obtained from a direct quantum 
mechanical calculation. 



11. THE LANGEVIN EQUATION 

The thermalization of the mode of frequency fl is described 
by the Langevin equation 



, + n'^u + Tu = £,{t)/fi 



(1) 



where u is the coordinate of the vibration; damping F = l/tr, 
and tr is the relaxation time; /i is the reduced mass of the 
mode; ^{t) is a normally distributed random force which de- 
scribes the interaction of the mode with the thermal bath of 
temperature T and the autocorrelation function 

imm) - 2firkBT^{t - 1') 

(the dimensionless function if{t) is normalized as 
/_ if{t)dt — 1), where ks is the Boltzmann constant. 

At thermal equilibrium the averaged energy of thermal vi- 
brations is defined by the relation 



+ 00 



rJo 2 

fi{uj^ + n^)\H{Lu)\'^ F{uj)duj , 



(2) 



where H{uj) ~ [fJ-i^'^ — uj'^ + icjF)] is the transmission 
function and F{u!) is the Fourier transform of the autocorre- 
lation function of the random force, 

f{lu) = ^ y (mm) eM-^^t}dt 



ip{t) exp{—iLut}dt 



(3) 



Thus, E — K + P where the averaged kinetic energy is de- 
fined by the relation 

1 ri .„ 

K = lim — / —fiii dt 

T^OO T Jq 2 



2kBTT 



uP'T{uj)dhj 



/o (f72 - Cj2)2 + ^2r2 

and the averaged potential energy is given by 

P= lim - / -^iVL^u^dt 



(4) 



2kBTT 



+ 00 



n^T{Lu)dLU 



(172 - tj2)2 + ^2r2 ' 



(5) 



where J-{uj) is the Fourier transform of the dimensionless au- 
tocorrelation function (p{t): 

1 r+°° 

T{lo) = -z- I f{t) exp{~iLLjt}dt . 
27r J_oo 

For a delta-correlated random force (the case of the white 
noise), Lp{t) = 6{t) and ^{u!) = l/27r. The integrals (|4|i 
and (|5]i can be easily calculated by a contour integration. The 
energy E = K = kBT/2 atV. ^ and E = K + P = kBT, 
K ^ P = kBT/2 for frequency fl > 0. 

For an exponentially-correlated random force (the case of 
low frequency color noise), (/3(i) — ^Xexp — \Xt\and!F{uj) = 
A2/27r((jj2 -|- ^2), where A = 1/ic and tc is the correlation 
time of the random force. In this case, the integrals (JUi and 
Q can also be calculated by a contour integration. The av- 
eraged kinetic energy K — kBTfxi^, F, A)/2 and the aver- 
aged potential energy P = kBTfp{Q., F, A)/2 with Jk and 
fp defined by 



fKin,T,x) 



A2(r!2 + A2-F2) + AFr>2 

(^2 + A2)2 - F2A2 



(6) 

(7) 



In the case of an harmonic oscillator, the Langevin equa- 
tion with white noise describes thermal vibrations of harmonic 
modes in the classical approximation, where the mean energy 
obeys E = kpT. In the case of a quantum harmonic oscilla- 
tor H = hfl{B+B + i), where h is the Planck constant, B+ 
and B represent creation and annihilation operators, the mean 
energy of thermal vibrations is given by 



E{Q,T) 



hfi 



exp{hn/kBT) - 1 



-hn 

2 



(8) 



The heat capacity of the oscillator is defined by c{fl,T) = 
dE{n,T)/dT ^ kBFE{n,T), where the Einstein function 
behaves according to 



FE{n,T) 



e^p{nn/kBT) 



hn 



kBT J [exp(arj/fcsT) - 1]2 



For T — > cxD, the Einstein function FE{fl,T) -^ 1, while in 
the limit T -^ 0, we get Fsi^, T) -^ 0. Consequently at high 
temperatures, the specific heat behaves as c{n, T) « kB and 
at low temperatures, one has c{n, T) « 0. For this reason, in 
the low temperature regime defined by T < Te = Tiil/kB, 
vibrations of the quantum oscillator are partially frozen. Con- 
sequently, a classical description of thermal vibrations is valid 
exclusively in the high temperature regime T > Te, where 
the Einstein temperature Te is defined by FE{n,TE) — 
e/(e- 1)2 = 0.9206735. 

If we drop the energy of vacuum vibrations hil/2, the ther- 
malization of the quantum oscillator can be characterized by 
the function 



G{n,T) = [E{n,T) - nn/2]/kBT 
nn/kBT 



exp{nn/kBT) - 1 



In the limit T ^ 0, the thermahzation coefficient G(f7, T) -^ 
0, while at T = Tb the function G{n, Te) = l/(e - 1) = 
0.5819767, and in the limit T -> oo, we have G{Q., T) -^ 1. 
Hence for temperatures T > 0, vibrations with frequency Vt > 
^e{T) — kBT/Ti will be frozen. We can thus conclude that 
it is uncorrect to model thermal fluctuations of these modes 
using the Langevin equation with white noise. 

As we stated at the beginning of this paper, the partial ther- 
mahzation of high frequency vibrations and the total thermal- 
ization of low frequency modes can be realized if one uses a 
Langevin equation with color noise which consists of low fre- 
quency components of the white noise. But it is necessary to 
consider in this case the temperature dependence of the noise 
correlation functions. For an exponentially-correlated random 
noise, this dependence can be deduced from the relation 



G{SIe{T), T) « [Jk(S^e{T),T, A) + fp{nE{T), r, A)]/2. 

In the limit T < VtE{T) = ksT/h, using Eq. Q and Q the 
last equation can be expressed in the simple form 



It is clear that the introduction of the color noise doesn't al- 
low to reproduce the temperature dependence of the quantum 
oscillator exactly. This is mainly due to the well known in- 
adequacy of a classical description to describe zero vibrations 
(ground state fluctuations). The figure nevertheless shows that 
the use of the color noise yields a qualitatively correct be- 
haviour of the heat capacity, that is, in the limit T — > the 
specific heat c ^ 0, and for T — > oo one has c ^ ks (with 
the increase of the temperature, the correlation time tends to 
zero and the color noise becomes a white noise). Most impor- 
tantly, we notice that there is a good agreement between the 
exact quantum result and the modified Langevin description at 
low temperatures T < Te where quantum effects dominate. 

We have shown that by using the Langevin equation with 
color noise, one can obtain a qualitatively correct picture for 
the temperature evolution of the specific heat in the case of 
a harmonic oscillator. The presence of non-linearities in the 
system will be tackled in the next section. 



III. THE EFFICIENCY OF THE PROPOSED SCHEME IN 
THE PRESENCE OF NON-LINEARITIES 



A2 



(9) 



x^ + ikBT/ny e-l' 

The equation (|9]l yields the following linear temperature de 
pendence of the correlation coefficient : 



A== 1/tc ^kBT/hy/e^^ 



(10) 



It follows from Eq. (fTol i that the description of the partial 
thermahzation of high frequency vibrations with the use of 
the Langevin equation ^ becomes possible if we introduce 
a correlation time tc that is inversely proportional to the tem- 
perature T of the thermal bath, i.e. 



tc = n^/F^^/ksT . 



(11) 



The time correlated color noise can be inserted in the 
Langevin description if one replaces the Langevin equation 
dTJ by the system of two equations. 



it = — i7^u — Fit + ^(t)//i, 

i = {vit)~m)/tc, 



(12) 
(13) 



where ?y(i) stands for the white noise generator, normahzed 
according to 



{7j{t)r]{t')) ^ 2fiTkBTSit - t'), 



As is well-known, anharmonicities are always present in 
physical systems and the validity of the harmonic approxi- 
mation which consists in representing the building blocks of 
a condensed system by linear oscillators is usually restricted 
to very low energy regimes. The anharmonic intermolecular 
forces may result either from the non-linearity of the individ- 
ual oscillators or from their non-linear mutual interaction. The 
natural question arises whether the Langevin equation with 
color noise can be used in the presence of non-linearities in 
the system. As a first attempt to answer this question, let us 
consider a non-linear oscillator whose dimensionless Hamil- 
tonian is given by 



H^lii' 



' ' ■ J^"^ 



r 



(14) 
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and tc is the correlation time whose temperature dependence 
is determined by (fTTT l. 

Fig. [T] compares the specific heat of a quantum harmonic 
oscillator and a classical oscillator whose dynamics is de- 
scribed by Langevin equations with color noise ( fT2] i and ( fTsT l. 



FIG. 1 : Temperature dependence of the heat capacity of the quan- 
tum harmonic (curve 1) and anharmonic oscillator (curve 3), classi- 
cal harmonic (curve 2) and anharmonic oscillator (curve 4) obtained 
from the Langevin equation with color noise (Q is the frequency of 
the oscillator and the anharmonicity parameter is chosen as /? = 0.1). 



The corresponding Langevin equation with color noise can be 
written in the form 



i = {V{t)-m)/tc, 

{v{t)^{t'))=2TT6{t-t'), 



(15) 



where T - the dimensionless temperature, the friction coeffi- 
cient r — 0.01 and the correlation time tc = \J& — 2/T. 

The numerical integration of the set of equations of motion 
(fTsT i yields the mean energy E — {H) of the anharmonic os- 
cillator as a function of T. Then the specific heat is computed 
from c{T) = dE/dT. 

In order to check the accuracy of the modified Langevin 
equation, we equally obtained the specific heat of this oscilla- 
tor by computing numerically the exact eigenvalues En of the 
quartic Hamiltonian (fl4] i. The diagonalization of the hamil- 
tonian matrix was performed in the basis of the harmonic os- 
cillator. The obtained eigenvalues were then used to find the 
partition function and the specific heat from the well-known 
relations 



Z{T) 

F{T) 

c{T) 



E' 



-E,JT 



-Tln(Z), 
d^F 



(16) 



-T 



dT"^' 



We compare in Fig. [T]the result of the numerical simula- 
tion to that obtained from the quantum statistical calculation 
(curve 3 and 4) at /3 = 0.1. We notice that the non-linearity ef- 
fects are indistinguishable within the accuracy of the proposed 
method. 

Let us note that Kleinert's variational path integral method 
allows us to obtain a fully analytical expression for the spe- 
cific heat of this quantum oscillator [2]. The method aims at 
approximating the quantum partition function 



Z = 



/'z)xe-^M, 



(17) 



where 



s^ I dT{^i' + y[x] 



(18) 



is the Euclidean action, with a quadratic trial action S'o[x] 
whose variational parameters are fixed by minimizing the 
right hand sight of the Jensen-Peirels inequality. 



F<Fo + T(S~So 



(19) 



This first order perturbation theory allows to compute the ap- 
proximative quantum partition function of a iV-body system 
as an integration over the classical configurational space. 



Z = 



dxo 



(27r/r)^/^ 



,-H'(xo)/T 



(20) 



where T4^(xo) is the so-called centroid potential and Xq stands 
for the centroid path [4 jj] . 

The effective potential W{xo) obtained by Kleinert for the 
trivial case of the quartic oscillator |2] is given in the ap- 
pendix. We verified that for /3 = 0.1, the analytical expression 
of the specific heat obtained from the partition function ( l20l l 
reproduces the numerical result (curve 3) of Fig. [T]with an 
error that is imperceptible on this scale. 

In order to examine the effect of non-linearities present in 
mutual interactions between vibrational modes, let us consider 
now the case of two linear oscillators with a dimensionless 
Hamiltonian 



H = -{ul+ul+ cjful + i^lul) + -/3(mi 



Ul) 



(21) 



where the variable Ui describes the dynamics of the i-th os- 
cillator (i = 1, 2), Wi the frequency of the same mode, and 
the parameter /? sets the strength of the non-linear interaction 
between the two oscillators. For this two-body hamiltonian, 
the system of Langevin equations with color noise takes the 
form 



U\ = —UJlUl + /3{U2 — Ul) —Till 

(2 = (r/2(t)-e2(t))Ac, 



-6(t), 



(22) 



where rii{t) stands for the random function that generates 
white noise with normalization conditions 

{m{tH{t')) = 2TTS{t-t'), z = 1,2, {viit)r,2it')) - 0. 

(T - dimensionless temperature, the dissipation coefficient is 
chosen as F = 0.01, and the correlation time is given by tc = 

It is also possible for this two-body system to obtain the 
temperature dependence of the specific heat from a purely 
quantum mechanical calculation, that is, by diagonalizing the 
hamiltonian matrix corresponding to Eq. (|2T] ) in the product 
basis of two harmonic oscillators of frequencies wi and uj2- 
The obtained energy levels were then used in Eq. ( |T6t in or- 
der to obtain the specific heat. 

For the sake of distinctness, we chose uji — I and uj2 = 10. 
The result of the numerical integration of the equations of 
motion (|22] ) illustrated in Fig. |2] indeed shows that the pres- 
ence of non-linearities in the interaction between individual 
modes improves the precision of the proposed Langevin ap- 
proach with color noise. To be more precise, for an interac- 
tion strength of (3 — 1, the temperature dependence of c{T) 
agrees better than the non-interacting case (3 = with the 
result obtained from quantum statistical mechanics. 




FIG. 2: Temperature dependence of the heat capacity for the system 
of two coupled linear oscillators (with frequencies LJi = 1,UJ2 = 10) 
obtained from a direct quantum mechanical calculation (curves 1, 3) 
and the Langevin equation with color noise (curves 2, 4). Solid lines 
correspond to decoupled oscillators (/3 — 0) and dashed lines to a 
coupling P — 1. 



teristic times for the integration of the generahzed Langevin 
equations requires some care. In fact, since the formula (fTTT i 
was obtained in the limit F 3> fcsT/Ti, the numerical value of 
the relaxation time should be large enough. It is thus crucial 
that the fixed value of the relaxation time remains always big- 
ger than the correlation time of random forces in the consid- 
ered temperature regime. On the other hand, very large values 
of tr are not suitable since it would require exceedingly long 
integration times to drive the system to thermal equilibrium. 
From a practical point of view, the numerical value tj, = 1 ps 
is adequate since the result remains practically invariant under 
the increase of this characteristic time. 

The proposed scheme will be applied in the next section to 
a Hamiltonian carbon nanotube model. 



IV. COMPUTATION OF THE HEAT CAPACITY OF 
CARBON NANOTUBES 



Using once more Kleinert's variational path integral 
method, we derived an analytical expression for the partition 
function of the two-body Hamiltonian (l2Tl i. The calculation 
is summarized in the appendix and the analytical result agrees 
very well with the numerical result in the temperature regime 
T > 0.3. 

These two examples clearly shows that the use of the 
Langevin equation with color noise can be used as a reason- 
ably accurate tool to compute the specific heat of complex 
molecular systems that are very difficult to study within the 
framework of quantum statistical mechanics. 

We will now present the scheme for evaluating the temper- 
ature dependence of the specific heat for general molecular 
systems by using a color noise. Let us define the N dimen- 
sional vector X = {a;„}^^]^ which denotes the spatial coordi- 
nates of the individual atoms in the system. In terms of these 
coordinates, the generalized Langevin equation describing the 
dynamics of the system takes the form 

Mx == -dH/dx-TMx + E, (23) 

s = (e-s)Ae, (24) 

where M - the mass matrix of the atoms, H - the Hamil- 
tonian of the system, the dissipation coefficient is given by 
F = l/tr iU - the relaxation time), Q ~ {'']n}n=i is the 
vector of normally distributed random forces obeying the nor- 
malization conditions 



{Vn{h)mit2)) = 2M„kBTSmS{h ~ h), 



(25) 



and the temperature dependence of the correlation time tc is 
still determined by the relation ( fTTT i. The choice of charac- 



Carbon nanotubes have the peculiarity of behaving as 
quasi-one-dimensional systems. This characteristic allows us 
to compute the specific heat of these structures using quantum 
statistical tools and this is what we will exploit in this sec- 
tion in order to check the efficiency of our modified Langevin 
approach on a concrete molecular system. 

For the sake of simplicity we will limit ourselves to the case 
of a nanotube with index of chirality (to, m). The structure 
of a carbon nanotube (CNT) with chirality (to, to) (armchair 
structure) is shown schematically in Fig. [3] The nanotube is 
characterized by its radius R, the angle shift (f and the longitu- 
dinal step h. The system consists of parallel transversal layers 
of atoms. In each layer, the nanotube has 2m atoms, which 
form m elementary cells separated by the angular distance 
A0 = 271 /m, so that h defines alternating longitudinal dis- 
tances between the transverse layers. Each atom of the CNT 
can be characterized by three indices (n, /, fc), where (n, I) de- 
fines an elementary cell (n = 0, ±1, ...,/ = 1, 2, ... , to), and 
k is the atom number in the cell, fc = 0, 1 (see Fig.[3]l. 

The Hamiltonian of the lattice of carbon atoms shown in 
Fig-EJcan be written in the following general form. 



^-im 



n 1 = 1 



iMiu. 



2 
■n,l,0 



+ <u)+^'M} 



(26) 



where M is the mass of a carbon atom, M = 

12X 1.6603 -lO-^^kg^Un^j^A; = {Xn,l,k{i),yn,l,k{t),Zn,l,k{t)) 

is the radius-vector that defines the position of the carbon 
atom {n,l,k) in the moment t, and the term Vn.i = 

P{Un-l,l,l, U„_i^( + i,o, Unjfi, Un^l^l, U„_|_i_i_i_i, Un+l.lfi) 

denotes the total potential energy given by a sum of three 
different types of potentials. 



n-l,l+2,0 n+l,l+l,0 

n,l+l,l 
n,l+l,0 



n-l,l-l,l 



n+l,l-2,l 




FIG. 3: Schematic representation of a carbon nanotube with chirality (m, m) and numbering of atoms in the structure. Thick red lines mark 
valent bond couplings, thin red arcs marks valent angle couplings, thin dashed lines show the foundation of two pyramids which form dihedral 
angles along the valent bonds in the elementary cell (n, I). For this figure, we chose m — 6. 



'Pn,l — ^(u„J,o, U„j^i) + V{Un-l,l+lfi, Unjs) + l^(u„,i,l , U„+i^i.o) + U {Un-iJ^i, Un^lfi, U„j^i) 

+ U{Unjfi, U„j^i, U„+i^;^o) + C^(Un-l,i+l,0, Un,l,l,Un+ld,o) + W^(UnJ,l , U„^;^o, Un-lJS, ^n+ld-is) (27) 

+ W{Unj^i,Undfi, U„+i^i_i^i, U„_ij^i) + W {Un-1,L1, ^71,1,0, U„^;^i , U„+i j-i^) + W{Un^ifi, U„_i^i , U„_i_i+i^o, U„+i^;_o) 
+ W{Un^lfi, U„j^i, U„+i^;^07 U„_i^; + i^o) + W^(u„--i J+i,o, U„_;^i, U„^i^o, U„+i^;^o)- 

I 



The first three terms describe a change of the deforma- with coordinates Ui and U2, characterized by the potential 
tion energy due to a direct interaction between pairs of atoms V^(ui, U2). The next six terms describe the deformation en- 



ergy of the angle between the links U1U2 and U2U3, taken into 
account with the potential [/(ui, U2, U3). Finally, the next 
six terms describe the deformation energy associated with a 
change of the effective angle between the planes U1U2U3 and 
U2U3U4, characterized by the potential VF(ui, U2, U3, U4). 

In our numerical simulations, we employ the interaction po- 
tentials frequently used in modeling of the dynamics of poly- 
mer macromolecules |4, 5, 6], 



F(ui, U2) = L'{exp[-a(p - pa)] - l}^ 



(28) 



where p = |u2 — Ui\, D = 4.9632 eV is the energy of the 
valent coupling, and po — 1.418A is the static length of valent 
bond; 



f/(Ui,U2,U3) =e„(cOS(^-COS¥3o)^ 



(29) 



where 



COSip = (U3-U2,U1 -U2)(|U3 - U2 1 • |U2 -Ui|) \ 

and cos 00 = cos(27r/3) = —1/2. Finally, 

VF(U1,U2,U3,U4) = et [l- (vi,V2)(|vi| • |V2|)"^] , 

(30) 
where Vi = (u2 — Ui) x (U3 — U2) and V2 = (U3 — U2) x 
(u4 — U3). The model parameters such as a = 1.7889 A^^, 
Eu = 1.3143 eV, and et = 0.499 eV can be determined from 
the phonon frequency spectrum of a planar lattice of carbon 
atoms 01]. 

Equilibrium structure of the nanotube with index (m, m) 
can be characterized by three parameters : its radius R, the 
angle shift ip and the longitudinal step h. Equilibrium posi- 
tions of the atoms in the tube are given by the coordinates 
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,1,0 


— h{n 


-1): 




<. 


1,1 


= /i(n 


-1). 






y'n. 


,1,0 


= RcOs{(j)ri, 


ll 


vt 


1,1 


= i?cos(0„,, 


i + ^\ 


(31) 


7^' 

n, 


1,0 ■ 


— i?sin(0„^, 


ll 


<„ 


!4 


= i?sin((/)„^; 


+ <P), 


1 



with cylindrical angles 0„^i = [I — 1 + {n — l)/2]A(f> and the 
angular distance Acj) = 2Tr/m. In order to find the parameters 
R, ip and h, we need to solve the minimization problem 



-f(Un-l,i,l' Un-l,i + l,0' "n,i,0' ^n,l,l''^ri+l,l-l,l,U^+i^Lo) 

— > min , 

R,ifi,h 



where we 

,,0 







have introduced the notations u! 

i^n,i.t^yn,i.t^z?i.i,t) and i = 0,1. The resulting value 
of the energy is then used as the minimum value. For a 
nanotube of the (6,6) type, we find a radius R = 4.1782 A 
and a longitudinal step h — 1.2590 A while for a nanotube of 
the (12,12) type, one obtains R = 8.3230 A and /i = 1.2560 
A. 

If one wishes to study small amplitude vibrations, it is more 
convenient to switch to local cylindrical coordinates Un,i,k, 

Vn,i,k, Wn,i,k defined by 



Xn,l,k — ^n,l,k ^~ '^n,l,k, 

yn,i,k = ?;°,_fc-w„,i,fcsin(/.°.;,s,+w„,i,fccos0" ; ,,, (32) 



Zn,l,k + Vn,l,k COS (j)°^ik + Wn,l,k sin 0° ; j., 
'"-1 



^n,l,k 

with the angle i^„.;.o = [I — 1 + (n ~ 1)/2]A0 and cpn,i.i — 
4'n,i,o + f- In this coordinate system, the Hamiltonian of the 
carbon nanotube takes the form 



m _. 



(33) 



n 1=1 



where the six dimensional vector x„.; — {un.i,o, Vn,i,o, Wn,i.o, Un,i,i, Vn,i,i, Wn,i,i) describes in local coordinates the shifting of 
the atoms located in the cell n, I from their equilibrium position. 
The equations of motion for the Hamiltonian ( |33] | are given by 

— MiCnj = ~Fn,l — Pxi(x„j;X„j+i;X„+i_i;X„+2,i_i;X„+2j) + Px2{^n,l-l',^n,r,^n+l,l~l',^n+2,l~-2',^n+2,l-l) 

+ -Px3(x„_l,;;X„_ij+i;X„_i;X„+i^;_i;X„+i_i) + Px^{Xn-2,l + l',^n-2,l+2',^n-l,l+l',^n,l]^n,l+l) (34) 

+ -Px5(Xn-2,i;x„_2J+i;x„_i^;;x„j_i;X„^;), 

with the function P^. == g|-P(xi, X2, X3, X4, X5), i — 1,2, ..., 5. In the Hnear approximation, the same set of equations takes 
the form 



— MXn,l — BiX„,l + B2Xn+l,l + i?2Xri-l,; + B^y.n+2,1 + -B3X„_2,; + -B4X„j+i + i?4X„^;_i 

+ -B5X„+ij_i + Bl:x.n-i,i+i + -B6X„+2,(-i + -BgX„_2,i+i + BjXn+2,1-2 + BjXn-2,1+2, (35) 

where the matrix coefficients are defined as 

Bl = ^XiXi + -Px2X2 + -Px3X3 + -fx4X4 + ^XsXs, -^2 = ^XiX3 + Px3X5 , ^3 — ^XiXs, 



XiX2 ~r -fx4X5 ; -O5 — -fx2X3 I -'x3X4i -^6 ' -'X1X4 I -'X2X5 7 ^^7 ~~ -'X 



r 



with the matrix of partial derivatives 



Xi ,XJ 



dxidxj 



(0,0,0,0,0), i,j = 1,2,3,4,5, 



The solution of the linearized equations (l35T l can be found in 
terms of plane waves in the form 



x„; = Ae exp(iqn + il5(j) — iojt), 



(36) 



where A stands for the amplitude of the wave, e - the unit 
vector of the amplitude, q e [0, tt] - the dimensionless wave 
number and 5 = 2TTJ/m (j — 0, 1, ..., m — 1) - the dimen- 
sionless orbital moment of the phonon. By injecting the ex- 
pression (l36b into the linear equation system ( |35] ). we obtain 
the eigenvalue problem 



Muj^A 



[Bi 



Boe'^ 



B*e-"i 



B.e^"^ 



r 



B*^-2iq 



Ble 



-B^e"^-'^ + Ble~"i+'^ + Bee^'"-'^ + B;e-^"'+'^ + Bje 



2iq-2iS 



-iS 



g*^-2rq+22S^^ 



(37) 



Thus, the calculation of the dispersion curves of the carbon 
nanotube requires the computation of the eigenvalues of the 
hermitian matrix of dimension 6x6 ( l37T i at each value of 
the wave number < g < tt and the moment S — 2Trj/m 
(j = 0, 1, ...,171 — 1). The dispersion curves obtained in this 
way consists of 6m branches (see Fig. |4](b)). 

The computation of the eigenvalues (l37b yields not only all 
the dispersion curves uj{q) but also the spectral density p{llj), 
normalized according to J„ p{u!)dijj = 1. 

A simple method for obtaining the temperature dependence 
of the spectral density consists in making a simulation of the 
thermal vibrations of the carbon nanotube. To this aim, the 
system is first driven to thermal equilibrium using the usual 
Langevin equations with white noise 



Mx„ 



TM±n,l + Sr; 



(38) 



with the dissipation coefficient T — 1/tr, tj. - the relax- 
ation time of atoms (it is reasonable to take tr = O.lps) and 
■=■«,; = (Cn,i,ij •■•jCn.i.e) - the six dimensional vector corre- 
sponding to the normally distributed random noises, describ- 
ing the interaction of the particles located in the cell (n, I) 
with the thermal bath. The characteristic correlation functions 
of these noises can be written as 

(Cn.fc.^(<l)6n,^,J(i2)) = 2A/fcsr,5„„,4;% (5(ti - ia), (39) 

where kB is the Boltzmann constant and T - the temperature 
of the thermostat. After having set the initial coordinates dSTT i 
and velocities to zero, the numerical integration of the sys- 
tem of Langevin equations was performed over t — 20 1^- 
We then decoupled the thermalized system from the bath and 
calculated the spectral density p{uj) of the kinetic energy dis- 
tribution of the atoms by following the real time dynamics of 
the isolated system. In order to increase the accuracy of the 
result, the spectral density was obtained from 100 indepen- 
dent thermalization processes and averaged over the atoms of 
the system. 



The spectral density profile obtained at T = 300 K is shown 
in Fig. |4](a). We notice that at this temperature, the density 
profile is in good agreement with the shape of the dispersion 
curves, which can be seen as a weak manifestation of non- 
linearity effects. If one assumes that the proper vibrational 
modes of the carbon nanotube remain linear, then the specific 
heat of the nanotube can be deduced from the integral 



c{T) 



Cq{uj) 



Cq{uj)p{uj)duj, 

exp(hLu/kBT) 
ksTj [exp{hu;/kBT) - 1]^ 



(40) 



huj 



where Cq{Lu) is the dimensionless thermal capacity of phonons 
with angular frequency uj. This method to compute the heat 
capacity of carbon nanotubes was first used in Ref. [?, UJ]. 

This approach is remarkably practical since the specific 
heat of the system follows from the simple knowledge of the 
spectral density p(w) of thermal vibrations . On the other 
hand, the spectral density can be deduced from the shape of 
the spectral curves, that is, by considering the later as the fre- 
quency spectrum of the harmonic modes in the nanotube. The 
heat capacity c{T) of the carbon nanotube with index (6,6), 
(12,12) and (10,0) obtained in this way is shown in Fig. |5] It 
is clearly seen in this figure that the specific heat is practically 
insensitive to the index of the nanotube. Furthermore, the tem- 
perature dependence of the specific heat is linear in the regime 
< T < 400 K. This temperature behaviour of the specific 
heat is in concordance with previous theoretical works 1191 UOll 
and confirms the experimental measures of the specific heat 
of titanium dioxide nanotubes 1 1 1 ] . 

We can equally obtain the spectral density that appears in 
Eq. (|40] | directly from dynamical simulations, i.e. by fol- 
lowing thermal vibrations of atoms in the carbon nanotube at 
finite temperature T > 0. As it can be checked in Fig. |5] 
numerical simulations show that these two approaches yield 
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FIG. 4: (a) Spectral density of thermal oscillations of armchair (6, 6) 
carbon nanotube for T = 300 K (red field corresponds to system 
equations with color noise) and (b) sixty dispersion curves of the 
phonon modes. Gray color marks the frequency region u < knT/h. 



almost the same result (the spectral density has a weak tem- 
perature dependence). 

Let us notice that the numerical integration of the Langevin 
equations with white noise dSSl) that allows to obtain the di- 
mensionless specific heat of the nanotube (m, 0) from the 
equation c{T) = {d{H)/dT)/{6NmkB), where Nh is the 
length of the nanotube, shows that the heat capacity is prac- 
tically independent of the temperature, that is, c = 1 over 
< T < 400 K (this results from the well-known equiparti- 
tion theorem of classical statistical mechanics, which states 
that the mean energy of each degree of freedom is equal 
to ksT). But the situation drastically changes, if one re- 
places the white noise of the Langevin equation with a time- 
correlated color noise whose temperature dependence is given 
by the formula ( fTTT ). In this case, thermal vibrations of the 
nanotube are described by the system of Langevin equations 



Mx„ 



^^n,l 






rMx„,; + s„ 



(41) 
(42) 



where 9„_; — (»7n,;,i, ■•-,??«, ;,6) - the six dimensional vec- 
tor corresponding to the normally distributed random noise 
and normalized according to ( [39] l, the relaxation time U = 




400 



FIG. 5: Temperature dependence of specific heat c of carbon nan- 
otube (6,6), (12,12) and (10,0). Blue line (curve 1) corresponds to 
the result obtained with the use of the frequency density of linear 
phonon waves while markers was obtained from the frequency den- 
sity of thermal vibrations. Red line (curve 2) corresponds to the result 
obtained from the numerical integration of Langevin equations with 
color noise. 



1 ps and the correlation time of random noises. 



The numerical integration of the equations of motion (|4TI) 
and (|42] | first yields the mean energy of the nanotube E = (H) 
versus temperature T. Then the specific heat is deduced from 
the relation c{T) = dE/dT. The result is illustrated in Fig. |5] 
First of all, it is clearly seen that the specific heat of nanotubes 
(6,6), (12,12) and (10,0) have practically the same tempera- 
ture behaviour, that is, it tends to zero for T -^ and it rises 
regularly with the increase of the temperature. Furthermore, 
the specific heat calculated with color noise coincides very 
well with the one obtained from Eq. (|40] | via computation of 
the spectral density. 

We show in Fig. |4l-(a) the shape of the spectral density 
of thermal vibrations of the carbon nanotube, obtained with 
color noise. We can note that only low frequency vibrations 
characterized by oj < ksT /fi are totally thermalized while 
the thermalization of high frequency modes is partial. More- 
over, the degree of thermalization decreases with increasing 
temperature. 

By considering the concrete example of carbon nanotubes, 
we have shown that the use of the Langevin equations with 
color noise W\\ and (l42T i allows to realize the quantum ef- 
fect of partial thermalization of high frequency modes and 
yields the correct temperature behaviour of the specific heat 
for complex molecular systems. The proposed method be- 
comes very useful especially for molecular systems with com- 
plex configurational dynamics, for example in the case of 
macromolecules which possess globular structure, or simply 
when one deals with a highly non-linear dynamics and it be- 
comes meaningless to consider the existence of a spectral den- 
sity of small-amplitude linear vibrations. 

Let us mention a further application of the generalized 
Langevin approach. The formalism that we have presented 
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consists in coupling to each particle a color noise of the same 
amplitude, which drives the systems to thermal equilibrium. 
If one instead couples the noise to chosen particles of the sys- 
tem or applies a noise of different amplitude to each particle, 
then the dynamics will be a non-equilibrium process. In this 
case, the use of the color noise may lead to new effects such as 
ratchet dynamics (see for ex. Ill2lll3ll ). which can be obtained 
neither within the usual Langevin approach with white noise 
nor with a quantum description of the dynamics. 



V. CONCLUSIONS 

We proposed a new method for computing the tempera- 
ture dependence of the heat capacity in complex molecu- 
lar systems. The proposed scheme is based on the use of 
the Langevin equation with low frequency color noise. We 
showed that the thermal behaviour of the correlation time of 
random forces, which is the key characteristic of the partial 
thermalization effect, can be described by a linear function 
of the inverse bath temperature, that is tc = h^/e — 2/kBT. 
We next illustrated non-linearity effects by considering two 
simple Hamiltonian models and we explicitly showed that the 
generalized Langevin approach can be used in the presence of 
anharmonicities in the Hamiltonian. Finally, by applying the 
proposed procedure to carbon nanotubes, we showed that the 
consideration of the color noise in the Langevin equation al- 
lows to accurately reproduce the temperature evolution of the 
specific heat in many-body systems. 

It is well-known that for complex systems having strong 
non-linearity effects in the quantum regime (T < Te), there 
exists a temperature gap unreachable by existing approxima- 
tive approaches such as the self-consistent phonon theory or 
the spectral density equations (|40] |. while the drawback of 
quantum Monte-Carlo methods is the large numerical cost in 
the case of many particle models. The proposed method may 
be very useful to fill this gap, especially if one wishes to inves- 
tigate the thermodynamics of realistic molecular systems with 
complex configurational dynamics, for example in the case 
of macromolecules which possess globular structures with a 
highly non-linear dynamics. 
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the first order variational path integral method. Since this ap- 
proach has already been intensively discussed [BbI] and ap- 
plied to several many body problems 11141 [Tsl UM . we will 
omit the technical details of the procedure and only report the 
analytical form of the effective potential for each Hamiltonian. 
The effective potential obtained by Feynman and Kleinert 
PO for the quartic quantum oscillator ( fT4l i can be written as 



W{x) 



?x^ + i(l + 3/3aV + ^«' 



3/3 
T 



^2 



Tin 



2T . f VL 

— — smn — 

n \2T 



, (Al) 



where the smearing parameter a and the frequency ft are de- 
duced from the following self-consistent equations : 



a = — coth — -, 

2n \2tJ n^' 

n'^ = 1 + 3(3(0^ +x^). 



(A2) 



These equations are solved in ref. fz] by a numerical iteration 
scheme at each x and T. Although this numerical procedure 
becomes very complicated when one deals with many body 
systems, it is shown in ref. [3] that the parameters a, H. and 
the centroid potential W{x) can be obtained for a one dimen- 
sional oscillator chain by Taylor expanding the equations ( IA2l i 
with respect to the non-linearity parameter /3 at the order 0(/3), 
then substituting the expansions in Eq. (lAll i and keeping only 
the terms of the same linear order, which finally yields the 
centroid potential in a fully analytical form. In the case of the 
quartic one-body potential (fT4l i. this procedure yields 



Wix) 



^x^ + f,{T)x^ + f,{T), (A3) 



where we have defined 



/i(T) ^ - -yr + — coth 



/2(T) 



3/3 
16 

Tin 



2T - coth 



2Tsinh 



1 

2T 
1 
2T 



1 

2r 

2 



(A4) 



The partition function that follows from Eq. ( IA3I ) and ( l20l ) 
can now be expressed in terms of the Bessel function of the 
second kind as 



APPENDIX A: THE SPECIFIC HEAT OF THE 
OSCILLATOR SYSTEMS FROM KLEINERT'S 
VARIATIONAL PATH INTEGRAL APPROACH 

We will give in this appendix the heat capacities of the sim- 
ple non-linear oscillator systems (fT4] i and ( 1211 1. obtained from 




(A5) 



and the specific heat is deduced from Eq. dTSI l. 

The effective potential of the two-body Hamiltonian (|2TI ) 
is computed in a similar way as for one dimensional oscillator 



chains yfl. The trial action that appears in Eq. ( [19] ) was chosen 
in the form 



dr i-{x^ 



y ) + Y^^' 



xo) 



+ -fiy-yo?+9{x~xo)iy-yo)\.(A6) 

where A^,, \y and g denote the trial parameters that minimize 
the Jensen-Peirels inequality ( fT9] l. After going through the 
usual minimization procedure and expanding the smearing pa- 
rameter a and the centroid potential at the order 0(/3), one ob- 
tains 



= - T 



X ' y 
^x y 



+ 



\ ujx coth I —^ I + ujy coth I —^ I > 

,w„ I \2T) ^ \2T) \ 



W{x,y) 



2uJx l~Oy 

Tin 



4T2 sinh(w^/2r) sinh(wy/2r) 



2+2^ + 4" 



+ '-{x-yf^%al{x-y)\ 



f' 
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(A7) 



The integration in Eq. (|20b can be carried out analytically by 
performing the coordinate transformation 



x-y 



x + y 



We finally obtain for the partition function 



(A8) 



Z 



' ' UJx'-'^y 



1 



X exp 



8y^7r/3T(L^2+^2)sinh(c^,/2r) sinh(^j,/2T) 
6/32a4 + ^4 



8/3T 



Ki 



8^ 



(A9) 



where we have defined 



f2 



2 _ '^^ ^y 



w:; 



- w,. 



3/3ag. 



(AlO) 
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